# Escoamento completamente desenvolvido em condutos de secao circular
#
# regime laminar (mu constante) ou turbulento (mu = mu(r))
# campo gravitacional nulo
# fluido newtoniano
# discretizacao por diferencas finitas
# 	EQUACAO DIFERENCIAL DE PRIMEIRA ORDEM
# diferencas centrais em pontos fora do contorno
# diferencas progressivas com tres pontos nos contornos com fdcoeffF.m
# code: FDM1pipeflow.m - octave
# para rodar: FDM1pipeflow OU run FDM1pipeflow.m

# Diego S. Rodrigues & Galera
# 23 de agosto de 2012
# PPG-CCMC - ICMC-USP
# Disciplina: SME-5802 - Introducao a Mecanica dos Fluidos Computacional
# Prof. Gustavo C. Buscaglia

clear all;

mu_lam = 1e-3;  # viscosidade dinamica do fluido
rho = 998.0;    # massa especifica
p0 = 1.0;  # pressao em z=0
pL = 0.0;  # pressao em z=L
L = 1.0;   # comprimento do tubo
r = linspace(0,1,101); # r(1) = 0.0, sempre!
gradp = (pL - p0)/L; # dp/dz
nunk = length(r) - 1; # numero de pontos (ou numero de incognitas - nunk)
R = r(length(r));
ystar = R/1000.0;
taus = gradp*R/2.0;
tol = 1e-6; # tolerancia para testar convergencia do metodo de Newton

for iunk=1:nunk
h(iunk) = r(iunk+1) - r(iunk);
end

for i=1:nunk
y(i) = R - r(i);
l(i) = R*(0.14 - 0.08*(r(i)/R)^2 - 0.06*(r(i)/R)^4);
if y(i) < ystar
l(i) = l(i)*(1-exp(-y(i)*sqrt(taus*rho))/(26.0*mu_lam));
end
end

u = ones(nunk,1);
uant = u;
erro = 1.0;

while (erro > tol)

# Construindo matriz jacobiana para metodo de Newton
for iunk=1:nunk
Res(iunk,1) = - mu_lam*u(iunk)/h(iunk) - rho*l(iunk)^2*abs(u(iunk))*u(iunk)/h(iunk)^2 - gradp*r(iunk)/2.0;
J(iunk,iunk) = mu_lam/h(iunk) + (rho*l(iunk)^2)/h(iunk)^2*(sign(u(iunk))*u(iunk) + abs(u(iunk)));
end

#Res = Res'; # transformando Res num vetor coluna

# Resolvendo o sistema linear
fprintf(stdout, "\n >>> Solving Linear System ... \n \n");
du = J\Res;
fprintf(stdout, "\n >>> Linear System was solved ... \n \n");

uant = u;
u = u + du;

erro = max(abs((u-uant)/uant)); # erro relativo usando a norma do maximo
printf("\n Erro = %e \n",erro);

end # end while

w(nunk+1) = 0.0;
for iunk=1:nunk
w(length(r)-iunk) = w(length(r)-iunk+1) - u(length(r)-iunk);
end

# Imprimindo saida num arquivo
fid = fopen("FDM1turb.dat","w");
fprintf(fid,"# r velw");

for i=1:(nunk+1)
fprintf(fid,"\n %g %g",r(i),w(i));
fflush(fid);
end

fclose(fid);

fprintf(stdout, "\n >>> END COMPUTATION \n \n");
